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The thermodynamics and dynamics of supercooled liquids correlate with their elasticity. In par¬ 
ticular for covalent networks, the jump of specific heat is small and the liquid is strong near the 
threshold valence where the network acquires rigidity. By contrast, the jump of specific heat and 
the fragility are large away from this threshold valence. In a previous work [Proc. Natl. Acad. Sci. 
U.S.A., 110, 6307 (2013)], we could explain these behaviors by introducing a model of supercooled 
liquids in which local rearrangements interact via elasticity. However, in that model the disorder 
characterizing elasticity was frozen, whereas it is itself a dynamic variable in supercooled liquids. 
Here we study numerically and theoretically adaptive elastic network models where polydisperse 
springs can move on a lattice, thus allowing for the geometry of the elastic network to fluctuate 
and evolve with temperature. We show numerically that our previous results on the relationship 
between structure and thermodynamics hold in these models. We Introduce an approximation where 
redundant constraints (highly coordinated regions where the frustration is large) are treated as an 
ideal gas, leading to analytical predictions that are accurate in the range of parameters relevant 
for real materials. Overall, these results lead to a description of supercooled liquids, in which the 
distance to the rigidity transition controls the number of directions in phase space that cost energy 
and the specific heat. 


I. INTRODUCTION 

Liquids undergo a glass transition toward an amor¬ 
phous solid state when cooled rapidly enough to avoid 
crystallization [1]. The glass lacks structural order: it 
is a liquid “frozen” in a local minimum in the energy 
landscape, due to the slowing down of relaxation pro¬ 
cesses. It is very plausible that the thermodynamics and 
the dynamics in supercooled liquids strongly depend on 
the microscopic structure of these configurations - here¬ 
after referred to as “inherent structures” [2]. However, a 
majority of glass theories [3-9] have focused on explaining 
the correlations between macroscopic observables seen in 
experiments (such as the relationship between thermody¬ 
namics and dynamics [10, 11]), while only a few [12-15] 
have investigated the role of structure. 

Experiments reveal that elasticity plays a key role in 
both the thermodynamic and dynamical properties in su¬ 
percooled liquids, such as the jump of specific heat and 
the fragility characterizing the glass transition. Specifi¬ 
cally, it has been found that (I) glasses present an excess 
of low-frequency vibrational modes with respect to Debye 
modes. The number of these excess anomalous modes, 
quantified as the intensity of the boson peak [16], shows a 
strong anti-correlation with the fragility [17, 18]. (II) The 
rigidity of the inherent structures is tunable by changing 
the fraction of components with different valences in net¬ 
work glasses [19-21], where atoms interact via covalent 
bonds and much weaker Van der Waals force. The co¬ 
valent network becomes rigid [22-24], when the average 
valence r exceeds a threshold Tc, determined by the bal¬ 
ance between the number of covalent constraints and the 
degrees of freedom of the system. Both the fragility and 
the jump of specific heat depend non-monotonically on 


r, and their minima coincide with rc [19, 25]. Interest¬ 
ing works using density functional theory [12, 26] inves¬ 
tigated the relationship between structure and fragility, 
but they do not capture this non-monotonicity. 

Recent observations [27-31] and theory [14, 32-40] in¬ 
dicate that in various amorphous materials, the pres¬ 
ence of soft elastic modes is regulated by the proxim¬ 
ity of the rigidity transition, linking evidence (I) and 
(II). To rationalize this connection, we have introduced 
a frozen elastic network model that bridges the gap be¬ 
tween network elasticity and geometry on one hand, elas¬ 
ticity and the thermodynamics and dynamics of liquids 
on the other [41]. This model incorporated the follow¬ 
ing aspects of supercooled liquids: (i) particles inter¬ 
act with each other with interactions that can greatly 
differ in strength, such as the covalent bonds and the 
much weaker Van der Waals interaction found in net¬ 
work glasses, (ii) Neighboring particles can organize into 
a few distinct local configurations, {in) The choices of 
local configurations are coupled at different location in 
space via elasticity. These features were modeled using 
a random elastic network whose topology was frozen, as 
illustrated in Fig. 2. The possibility for local configura¬ 
tions to change was incorporated by letting each spring 
switch between two possible rest lengths. Despite its sim¬ 
plicity, this model recovered (I) and (II). In particular, 
it reproduced the non-monotonic variance of the jump of 
specific heat and the fragility with the coordination z of 
the network: they are extremal at Zc = 2d {d is the spa¬ 
tial dimension), where a rigidity transition occurs. This 
model could be solved analytically, and it led to the view 
that near the rigidity transition, the jump of specific heat 
is small because frustration vanishes: most directions in 
phase space do not cost energy, and thus do not con- 
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tribute to the specific heat. 


Rigid 



floppy isostatic stressed 


FIG. 1: (Color online) Illustration of rigidity transition. 
Blue, green, and red color the floppy, isostatic, and stressed 
clusters, respectively. 


thermodynamic predictions of Ref. [41] relating struc¬ 
ture to the jump of specific heat. In this work, we di¬ 
rectly show numerically and theoretically that the predic¬ 
tion for the jump of specific heat is essentially identical 
in adaptive and frozen elastic network models. Section 
II describes the adaptive network models. Section III 
presents the numerical results of the model, while Sec¬ 
tion IV gives the explicit derivation of the thermody¬ 
namic properties, developing an approximation scheme 
to deal with the temperature-dependence of the number 
of over-constraints in the system, treating them as an 
ideal gas. 


This is a novel explanation for a long-standing prob¬ 
lem, and it is important to confirm that this view is ro¬ 
bust when more realism is brought into the model. In 
particular, the model used frozen disorder to describe 
elasticity, whereas it is itself a dynamical property in liq¬ 
uids, where there cannot be any frozen disorder. The 
thermal evolution of the topology of the contact net¬ 
work and its effects on rigidity transition were also not 
addressed. A network is rigid when an imposed global 
strain induces stress, and the rigidity can be achieved 
topologically by adding constraints [22], see Fig. 1 for an 
illustration in a small network. The network is said to 
be self-stressed if some of the constraints are redundant, 
removing those leaves the network rigid. Three scenarios 
of rigidity transition have been extensively studied in the 
literature [42, 43] (but see Ref. [44] for a recent fourth 
proposition). Spatial fluctuations of coordination are im¬ 
portant in the first two. The rigidity percolation model 
[45-48] assumes that bonds are randomly deposited on 
a lattice. Fluctuations lead to over-constrained (self- 
stressed) clusters even when the average coordination 
number is not sufficient to make the whole network rigid. 
This model corresponds to the infinite temperature limit. 
To include these effects, self-organized network models 
were introduced [49-52], where overconstrained regions 
are penalized. A surprising outcome of these models is 
the emergence of a rigidity window: rigidity emerges at a 
small coordination number before the self-stress appears 
(even in the thermodynamic limit). Finally, in the mean- 
field or jamming scenario, fluctuations of coordinations 
are limited. Similar to the simple picture in Fig. 1, the 
rigidity, and the stress appear at the same Zc in the ther¬ 
modynamic limit. The rigid cluster at Zc is not fractal 
and is similar to that of packings of repulsive particles. 
The model of Ref. [41] assumed that networks were of 
this last type. 

Recently, we have introduced adaptive elastic network 
models [42], where the topology of the network is free to 
evolve to lower its elastic energy as the system is cooled. 
We found that as soon as weak interactions are present, 
the network of strong interactions becomes mean-field 
like at low temperature. However, the thermodynamic 
properties were not studied to test the robustness of the 


II. MODEL 



FIG. 2: (Color online) (a) and (b) Illustration of the frozen 
network model [41]; (c) and (d) illustrate the adaptive 
network model [42]. In the latter case, the triangnlar lattice 
is systematically distorted in a unit cell of four nodes shown 
in the inset of (c). We group nodes by four, labeled as. A, B, 
C, and D in Fig. 2. One group forms the unit cell of the 
crystalline lattice. Each cell is distorted identically in the 
following way: node A stays, while nodes B, C, and D move 
by a distance <5, B along the direction perpendicular to BC, 
C along the direction perpendicular to CD, and D along the 
direction perpendicular to DB. 5 is set to 0.2 with the lattice 
constant as unity. Weak springs connecting (b) six nearest 
neighbors without strong springs and (d) six 
next-nearest-neighbors are indicated in straight cyan lines, 
emphasized for the central node, (c) Illustration of an 
allowed step, where the strong spring in red relocates to a 
vacant edge indicated by a dashed blue line. 
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In our model degrees of freedom are springs, which are 
poly-disperse and can move on a lattice. The lattice is 
built using a triangular lattice with periodic boundary 
conditions, see Fig. 2(c), with a slight regular distortion 
to minimize the non-generic presence of zero modes that 
occurs when straight lines are present, as illustrated in 
the inset of Fig. 2(c). Polydisperse and mobile “strong” 
springs of identical stiffness k connect the nearest neigh¬ 
bors on the lattice and model the covalent constraints. 
We model weak Van der Waals interactions with “weak” 
and stationary springs of stiffness few “C fc adding to 
all next-nearest-neighbors on the triangular lattice, il¬ 
lustrated in Fig. 2(b). We introduce a control parameter 
a = {z^/d){k^/k) to characterize the relative strength 
of the weak interactions, where the spatial dimension is 
d = 2 and the number of weak constraints per node is 
chosen Zw = 6 . 

The number of “covalent” springs W, equivalent to 
the coordination number z = 2 Ns/N {N is the number 
of nodes in the lattice), is also a dimensionless control pa¬ 
rameter. For a given 5z = z — Zc, the valid configurations 
are defined by the locations of the Ng springs, indicated 
as F = {7 o where the Greek index 7 labels 

springs and the Roman indices (*,j) label the edges on 
triangular lattice between nodes i and j. We introduce 
the occupation of an edge: cr^j = 0 if there is no strong 
spring on the edge ij, and = 1 if there is one. If 

denotes the geometric length between nodes i and 
j on the lattice, we assume that the spring 7 has a rest 
length Ij = r(ij) +£ 7 , where the mismatch is a feature 
of a given spring, are sampled independently from a 
Gaussian distribution with mean zero and variance e^, 
which thus characterizes the polydispersity of the model. 
ke^ is set to unity as the natural energy scale. 



FIG. 3: (Color online) Illustration of configuration energy 
of the adaptive network model {Sz = 0.27). Solid lines are 
springs, colored according to their extensions: from red to 
purple, the springs go from being stretched to being 
compressed, with spring extensions shown in the unit of e. 
Left: Nodes sit at lattice sites, so the color shows the rest 
length mismatches of the springs {£7}. Right: Nodes are 
relaxed to mechanical equilibrium. Most links appear in 
green, indicating that most of the elastic energy is released. 
The configuration energy is dehned by the residual energy. 


The energy of an inherent structure is denoted 'H(r). 
The configuration F is sampled with probability propor¬ 
tional to exp(—'H(F)/T) in the liquid phase, with ks = 1. 
Temperature T serves as a third dimensionless control 
parameter. H(F) is defined as the remaining energy once 
the nodes of the network are allowed to relax to mechan¬ 
ical equilibrium: 


n{T) 



\\Ri Rj 



<hl>2 



( 1 ) 


where Ri is the position of particle i and (*, j )2 labels the 
next-nearest neighbors. The minimal energy can be cal¬ 
culated by steepest decent as illustrated in Fig. 3, but this 
is computationally expensive. Instead, we approximate 
the elastic energy in the linear response range, setting 
that <C 1 [53]. The above minimization expression 
Eq.(l) could then be written as, 

'^(r) = ^ {i,i) ,{l,m)^{l,m) + o(c ) (2) 

r 

where = £7 when spring 7 connects i and j. The 
coupling matrix Q — 'P — S{S*S + derived 

in our previous works [41, 42] (or see Appendix Sec. A), 
is a product of the structure matrix S and its transpose 
S*, the structure matrix of the weak spring network iSw, 
and V the projection operator of the triangular lattice 
onto occupied edges. The structure matrices S and 5w 
describe the topology of the networks of strong and weak 
springs: if neighbor nodes i and j are connected, the 
change of the distance between i and j, ■ 

SRi+S^ij)j-SRj+o{SR^), due to displacements of nodes 
SR. We point out that as the weak network is fixed, S 
and thus G depend only on the network topology of strong 
springs, but not on the mismatches £7. 


Our model is a generalization of on-lattice network 
models: setting the interaction strength control param¬ 
eter a = 0 , it naturally recovers the randomly diluted 
lattice model [48] when T = 00. It is also related to 
the self-organized lattice model [49, 50], which postu¬ 
lates that elastic energy is linearly proportional to the 
number of redundant constraints [49, 54]. We will find 
that this assumption holds true for a = 0 and T <C 1 . 
However, the existence of weak interactions among sites 
means that in real physical systems o; > 0. This turns 
out to completely change the physics, an effect that our 
model can incorporate. 
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III. NUMERICS 

We implement a Monte Carlo simulation to sample the 
configuration space of the model, with 10® Monte Carlo 
steps at each T. At each step, a potential configuration 
is generated by a Glauber dynamics - moving one ran¬ 
domly chosen spring to a vacant edge, as illustrated in 
Fig. 2(c). We numerically compute the elastic energy of 
the proposed configuration using Eq.(2): calculating the 
structure matrix S and then the corresponding G- On 
computing G, the matrix inversion, (5*5 -I- ^5^5w)“^, 
is singular when the network contains floppy structures, 
which do not appear except when few = 0. When a = 0, 
we implement the “pebble game” algorithm [55] to iden¬ 
tify the over-constrained sub-networks, and then do ma¬ 
trix division in the subspace, as the isostatic and floppy 
regions store no elastic energy after relaxation. We have 
found little finite size effect by varying the system size 
from iV = 64 to iV = 1024 nodes in the triangular lat¬ 
tice. In the following, we present our numerical results 
of networks with N = 256 nodes, averaged over 50 real¬ 
izations of random mismatches if not specified. 


A. Dynamics 

We investigate the dynamics by computing the correla¬ 
tion function ^(t) = Af^(i_j^/^/ 3 jv) ((gWk(O)) - 
where |cr(t)) is the vector indicating the occupation of 
all edges at time t. The correlation C{t) decays from 
one to zero at long time scales. We define the relaxation 
time T as the time C'(t) = 1/2, and the numerical results 
of r as a function of temperature T for several different 
coordination numbers are shown in the Fig. 4. 



FIG. 4: (Color online) Relaxation time r in log-scale versus 
inverse temperature 1/T for different coordination numbers 
Sz and a = 0.0003. The solid black line indicates a power 
law relation between r and T: t ^ . 

We find that the implemented dynamics is not glassy. 



FIG. 5: (Color online) Left: Shear modulus of adaptive 
networks at temperature T rescaled by G at T = oo 
G(2, T)/G(a, oo), a — 0.0003. The temperature T is rescaled 
by Tg. Right: Correlation between transition temperature 
Tg and shear modulus G in the frozen network model [41]. 


The relaxation time increases as a power law of the tem¬ 
perature even much slower than a strong glass 

that would display an Arrhenius behavior log^o t oc 1/T. 
This result is very surprising because the frozen elastic 
network model we studied earlier was glassy (its fragility 
was similar to that of network liquids). Despite being 
dynamically very different, these two models are almost 
identical as far as thermodynamics is concerned, as we 
will see below. It could be that the lack of glassiness 
comes from our choice of Monte-Carlo where springs can 
try other locations anywhere in the system [56]. 

To compare the thermodynamics of these models we 
now need to define an effective glass temperature Tg 
(even if we do not see a real glass transition). We 
do that by using the empirical Lindemann criterion [57] 
according to which an amorphous solid melts when the 
standard deviation of particles’ displacements is 

greater than a fraction cl of the particle size a. The co¬ 
efficient Cl must depends on the quench rate q, since this 
is also the case for Tg. This dependence is logarithmic, 
because the dependence of relaxation time on tempera¬ 
ture in experimental glass formers is at least exponential 
(for typical experimental quench rate in supercooled liq¬ 
uids, Cl ~ 0.15 [58]). We can estimate this standard 
deviation via the elastic modulus if we treat the glass 
as a continuum {SB?) ^ T/Ga where G is the instan¬ 
taneous shear modulus of the structure [8], we thus get 
Tg oc Ga^/\n{l/q). We set the lattice length a in our 
model to unity. 

We measure the shear modulus averaging over config¬ 
urations at given temperatures, shown in the left panel 
of Fig. 5. Practically, we choose Tg = (G)Tg/ln(l/10®g), 
where the cooling rate q is defined as the inverse of the 
number of Monte Carlo steps performed at each tempera¬ 
ture in the model. (•)Tg is the mean value at temperature 
Tg. The prefactor in this definition of Tg does not affect 
qualitatively our conclusions, but for this pre-factor the 
definition of Tg in the frozen model [41] is essentially iden- 
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tical to the dynamical definition used in [41], as shown 
in the right panel of Fig. 5 by lining up G and Tg. The 
specific values of Tg following that definition are shown 
in the inset of the bottom panel of Fig. 7, they corre¬ 
spond to Tg = (G)Tg/ In(lO^) in the present model, and 
Tg = (G)Tg / ln(10®) in the frozen network model [41], 
which is simpler to simulate and can thus be equilibrated 
longer. 


B. Specific heat 

The specific heat data shown in Figs. 6 and 7 are our 
central numerical results. The energy E = {%) is ob¬ 
tained using a time-average over Monte Carlo steps, and 
is shown in Fig. 6(a). The specific heat is calculated as its 
derivative c = -^dG/dT, and is shown versus T for sev¬ 
eral coordination numbers when o = 0 in Fig. 6(b) and 
a = 0.0003 in the top panel of Fig. 7. When a = 0, the 
specific heat increases as temperature decreases for net¬ 
works with Sz > 0 while it meets a maximum at Ta ~ 1 
and decreases under cooling when T < Ta ii 6z < 0. By 
contrast, the specific heat increases under cooling close to 
the transition temperature for all coordination numbers 
when a > 0. In addition, when T < a, c ^ 0.5.All these 
results are qualitatively identical to our previous frozen 
model. 

To define the jump of the specific heat at the glass 
transition, we simply measure the specific heat at our 
glass transition Tg defined above. This definition is nat¬ 
ural, since in a real glassy system, below Tg the liquid is 
essentially frozen in an inherent structure, and the contri¬ 
bution to the specific heat from configurational entropy 
(i.e. the bottom energy of inherent structures) vanishes. 

Our central numerical result is shown in the bottom 
panel of Fig. 7: c(Tg) varies non-monotonically with the 
coordination number z when a > 0. When the network 
of strong springs is poorly coordinated ;< 0, c(Tg) de¬ 
creases as z increases; When the strong network gets bet¬ 
ter coordinated Sz > 0, c gradually changes to increase 
with z; c is minimal at the proximity of the rigidity tran¬ 
sition for finite a. These numerical results are very 
similar to empirical observations, see Point (II) in the 
introduction. Our data are in fact very similar to that 
of the frozen model, which essentially follows the dotted 
lines in Fig. 7. 


C. Number of redundant constraints R 

When 0 = 0 and T —>• 0, the specihc heat is simply 
proportional to R, as shown in Fig. 6(b). This number 
is fixed, R = N6z/2, in the frozen network models. It 
varies in the adaptive network model and depends on the 
temperature. As the Maxwell counting gives the minimal 
number of redundant constraints of a network, we can 


define an excess number of redundant constraints 

nex^^(A-^0(<5z)), (3) 

where 0(a;) is the Heaviside step function. Uex counts 
the average number of redundant constraints, additional 
to the Maxwell counting. This excess number of redun¬ 
dant constraints decreases monotonically to zero under 
cooling. When a = 0, riex is proportional to Vt in the 
adaptive network model at low temperature, shown in 
Fig. 6(c). 


IV. THEORY 


As illustrated in Fig. 8, in the frozen elastic model 
we found that as a —?► 0, c converges to a constant if 
z < Zc, whereas it behaves as z — Zc for z > Zc. As a 
is increased, the discontinuous behavior becomes smooth 
and looks similar to experimental data. We seek to derive 
these same features in the adaptive network models. 


A. Thermodynamics 

For simplicity, we consider the annealed free energy 
.^ann = —TliiZ. It is exact in the random energy 
model [59] above the ideal glass transition [60] and we 
find it to be a good approximation of J- in our mod¬ 
els [41]. The over-line implies an average over disorder 


^ = E E exp[-7f(r)/T] (4) 

{a} perm[7] 


where a given configuration F is characterized by {a} 
indicating which edges are occupied on the triangular 
lattice, and perm[ 7 ] labels the possible permutations of 
springs’ rest lengths. 

We first average over the quenched randomnesses. Us¬ 
ing the linear approximation Eq.(2) and the Gaussian 


distribution p(e.y) = 






^ = Ef^)!exp 


{'^1 


- - tr In I I - 


GiW}) 

T 


( 5 ) 


The factorial comes from Ngl = X]perm[ 7 ] 1 as t/ is inde¬ 
pendent of the permutation. I is a 3N x 3 A identity ma¬ 
trix; each component corresponds to an edge on the lat¬ 
tice. To compute the trace in the exponent, we first make 
the approximation that the weak springs are weak and 
numerous ~ ^TjsidxNd, which corresponds to the 

highly connected limit Zw —t oo and finite a. We can then 
decompose the coupling matrix Q 'P—S{S*S+aI)~^S^ 
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FIG. 6: (Color online) Thermodynamics of the adaptive network model without weak constraints a = 0. (a) Energy E/Ns vs 
temperature; (b) Specific heat C/Ns vs temperature; (c) Excess number density of redundant constraints nex extracted using 
the pebble game algorithm vs temperature. Symbols are numerical data, solid lines are theoretic predictions. 


as [41]: 

IV'p)(V'p|+ (6) 

p({(t}) l^({cr})>0 

where p labels the vectors \tpp) satisfying = 0 


(i.e. a basis for the kernel of 5*), and where the \tpui) 
satisfy SS*\tpui) = The number of redundant 

directions is 1 = Ns — {Nd — F) = R. Note that 
tr?^ = Ns, Nd — F gives the number of frequencies w, 
and F counts the number of floppy modes. The modes 
l^p); R, and uj depend on occupation {cr}. As the 

\ip)^s are orthonormal, the trace in Eq.(5) gives 


Z = 



! ^ exp 


r,D{u) 


Ns 


(^s{nr,D{uj)) 


■ln(l+y) 



du}D{uj) ln(l 


1 a 

TuFTa 



(7) 


where s(nr,D{w)) = lnX]{a} is configura¬ 

tional entropy density with given number of redundant 
constraints tir = R/Ns and density of vibrational modes, 
D{lo), satisfies (1 - Ur) f dujD{uj) = limjv-i-oo Z]c^> 0 ' 

B. No weak interactions 

Neglecting the weak constraints a = 0, the last term in 
the exponential vanishes and the summation over states 
with given density of states can be absorbed into the 
entropy, which then depends only on the number of re¬ 
dundant constraints. 

Z = f in(i+^)] (g) 

'' rir 

We propose an ideal-gas picture of “defects” to find 
an approximation form of the entropy s{nr). When the 
coordination number is very small z < Zc and the net¬ 
work is mostly floppy, redundant constraints are defects 
localized in rigid islands. Similarly, when the coordina¬ 


tion number is very large z > Zc with most regions of 
the network rigid, there are localized floppy modes in 
regions where there are negative fluctuations of coordi¬ 
nation number, which we again described as defects, see 
illustration in Fig. 9. The number of such floppy modes 
is equal to the number of additional over-constrained in 
the rigid cluster. The entropy gains from having these 
defects. Assuming that such defects are independent, we 
approximate the entropy by that of an ideal gas: 

s(»^ex) ~ s(0) - nex In , (9) 

eno(2;) 

where nex is the excess number of redundant constraints 
defined in Eq.(3) and is thus counting the number of 
defects. s(0) is the entropy density of the states with 
a minimal number of redundant constraints (i.e. they 
satisfy the Maxwell counting); and no{z) is the excess 
number of redundant constraints at T = oo. Both s(0) 
and uq depend only on z and the lattice structure. This 
form of Eq.(9) fails when the assumption of independent 
“defects” breaks down, as must occur near the rigidity 
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FIG. 7: (Color online) Top: Specific heat c{z,T) vs scaled 
temperature TjTg for networks with average coordination 
numbers near and away from the isostatic on both floppy 
and rigid sides. The strength of the weak constraints is 
given by Q = 0.0003. Bottom: Specific heat at temperature 
Tg, c{z,Tg), vs coordination number Sz for 
a = 0, 0.0003, 0.003, 0.03. The inset shows the transition 
temperature Tg for different a and a. Symbols are numerical 
results, and lines are theoretical predictions: dashed lines 
are for frozen network model and solid lines are for the new 
model derived in section IV. 


transition. However, our numerical results indicate that 
this approximation is very accurate, we see deviations 
only for jJzl < O.I. 

We numerically test the formula Eq.(9) for a triangular 
lattice. The configurations with R redundant constraints 
are weighted by for different values of the parameter 
p. From Eq.(9), the mean and variance of the excess 
number density of redundant constraints, riex, satisfy the 
following formulas: 

d s 

P=-z - ^ nex(z,/3) = no(z)e“^ (10a) 

t/TT-ex 



FIG. 8: (Color online) Theoretical predictions for the jump 
of specific heat. For vanishingly weak springs a ^ 0, it is 
predicted that the jump is essentially constant for z < Zc 
and then drops to zero a Zc- For larger z, it behaves as 
2 — Zc- As a grows this sharp curve becomes smooth, but a 
minimum is still present near z = Zc- 



FIG. 9: (Color online) (a) 2 < 2c, localized redundant 
constraints (red) in a floppy sea (blue); (b) z > Zc localized 
floppy modes (blue) in a rigid sea (red and green). 



FIG. 10: (Color online) Left: Excess number density of 
redundant constraints nex(2,/3)- Right: Fluctuation of the 
number density of redundant constraints (Aric)^. The solid 
black lines show the predictions from the approximate 
entropy Eq.(9). 


P) 


-0^^n^^{z,P) = P‘^ne^iz,P) 


(10b) 


Our numerical results coincide with Eqs.(lOa) and (10b) 




























































































remarkably well, with minor deviations for \Sz\ ^ 0.1, as 
shown in Fig. 10. 


Applying Eq.(9), we derive the thermodynamics of our 
model when a = 0. Solving the saddle point of Eq.(8), 
we obtain the average energy density: 


1 , N fo + n'exiz,T) T 

2 ITT 


(11a) 


the specific heat: 

^-c(z,r) = 


^0 + ^rie^{z,T) 


1 


(11b) 


2 (1+T)2 

and the excess number density of redundant constraints: 


/ iX-l/2 

nex(z, T) = no(z) ( 1 + - j (11c) 


C. General case 


In the thermodynamic limit, Ng —^ oo, we take the 
saddle point of Eq.(7), 


|^ = ln(l + i)-/<la,DMln(l + 

and 

25s M 1 OL 

—— = (1 - n^) In 1 + - 
oDyuj) \ i a;2 + a 

and solve for energy, 


1 a 

Tuj^ + a 


(12a) 


(12b) 


N.. 


■E{z,T,a) 


nr{T) T 
2 l + T 


1-nJT) 


dujD{ui, T) 


oT 

a + (a;2 + a)T 


(13) 


where tq = ^■Q{5z). 

As no( 2 :) is expected to be an analytic function of z, 
Eqs.(ll) indicate that c converges to the one found in 
frozen network model in the limit T —>■ 0: c = 0 when 
Jz < 0 and c = (5z/2z when Sz > 0 - the dashed yellow 
line in Fig. 8. This is our first central result, which shows 
that our previous results hold even when the network is 
adaptive. 

Eqs.(ll) predict the energy, specific heat, and the num¬ 
ber density of redundant constraints at an arbitrary tem¬ 
perature without any fitting parameter. The solid lines, 
shown in Fig. 6(a) and (b), are predictions of Eqs.(lla) 
and (11b), respectively, with riex as the numerical in¬ 
put. They are closely consistent with the data points, 
which confirms the annealed free energy approximation 
when a = 0. A power-law with numerical prefactor 
no(z) = nex(z,oo) predicted by Eq.(llc) coincides well 
with data points in Fig. 6(c). 

Extending to finite glass transition Tg at a = 0, we find 
a correction vanishing as in addition to c ~ i5z/22, 
assuming ^ G ~ (5z for z > Zc- But this correction 
is quantitatively unimportant as uq < 0.03 and does not 
change qualitatively the linear growth of the specific heat 
when 6z > 0, as illustrated by the solid orange line in 
Fig. 8. 

Our theoretic prediction that riex 0 when T —^ 0 
validates the assumptions of [49, 50, 54] that the energy 
of redundant bonds is proportional to their number, and 
that this number is i?o at T = 0. 


The specific heat predictions from differentiating Eq.(13) 
with numerical inputs nr-(z, T, a) and Dz,T,a{uj) are plot¬ 
ted as solid lines in Fig. 7. (See Appendix Secs. BCD 
for the temperature dependence of D{uj).) Notice 
that replacing nr(T) by 5z/z and D{uj,T) by its low- 
temperature limit D{uj) studied in [41, 61, 62], Eq.(13) 
recovers exactly the one obtained in the frozen network 
model, whose predictions are plotted as dashed lines in 
Fig. 7. The dashed lines converge to the solid lines de¬ 
spite differences at high temperatures for weakly coordi¬ 
nated networks. 

In the limit a —>■ 0 and T <C a, Eq.(13) converges to 
E/Ns = T/2, which indicates a constant specific heat 
c = 0.5 when 5z < 0 independent of the models. This is 
shown by the solid orange line and the dashed yellow line 
in Fig. 8, and is our second key theoretical result showing 
the robustness of our conclusions for adaptive networks. 


V. CONCLUSIONS 

In this work, we have studied the correlation between 
the elasticity of inherent structures and the thermody¬ 
namics in covalent glass-forming liquids using adaptive 
network models. We found numerically and explained 
theoretically why these models have a thermodynamic 
behavior similar to frozen network models [41] which cap¬ 
tures nicely experimental facts. 

The main prediction conclusion of [41] is thus robust: 
as the coordination number approaches Zc from above, 
elastic frustration vanishes. This leads both to an abun¬ 
dance of soft elastic modes, as well as a diminution of 
the number of directions in phase space that cost energy, 
which is directly proportional to the jump of specific heat. 
Below the rigidity transition, the elasticity of strong force 
network vanishes, thus the energy landscape is governed 
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by the weak Van der Waals interactions. At these energy 
scale, all directions in contact space have a cost, and thus 
the specific heat increases. Thus thermodynamic prop¬ 
erties are governed by a critical point at i52: = 0 , a = 0 
where the jump of specific heat is zero. This prediction 
focuses on the configurational part of the jump of specific 
heat, since we considered only the energy minima in the 
metastable states. In Appendix Sec. E, we argue that the 
vibrational contribution to this jump is so small in our 
models. Thus the main prediction of the specific heat 
still holds, even when including the vibrational part. 

Beyond network glasses, our main result potentially 
explains the correlation between elasticity and the key 
aspects of the energy landscape in molecular glasses [19, 
25, 63]. Indeed according to our work we expect glasses 
with a strong Boson peak to display less elastic frustra¬ 
tion, so that they have a limited number of directions in 
phase space costing energy, see discussion in [41]. 

We thank E. DeGiuli, G. Diiring, J. Lin, E. Lerner, 
G. Sandford for discussions, and D. Jacobs for sharing 
the pebble game code. This work has been supported 
primarily by the National Science Foundation Grant No. 
GBET-1236378, and partially by the Sloan Fellowship, 
the NSF Grant No. DMR-1105387, and the Petroleum 
Research Fund Grant No. 52031-DNI9. 


as indicated in Eq.(Al). Minimizing Eq.(l), one gets: 


H{T) 


min 

{SR.} 


2 S^.i5Ri + 

7 i 


-t- 


= min 


<7 i ) 

^ \{e\V\e) + 2{e\S\6R) + {SR\M\SR) 


(A2) 


where we use bra-ket notations to indicate summation 
over edges or nodes, V projects the edge space to the 
subspace occupied by springs, M = S^S + is 

the stiff matrix connecting the responding forces and dis¬ 
placements of nodes in an elastic network [64], and •* 
is our notation for the transpose of a matrix. Solving 
Eq.(A2), one finds the linear response, 

\SR) =-M-^S*\e) (A3) 


which for a given mismatch field je) minimizes the elastic 
energy in Eq.(l). Inserting Eq.(A3) back into the linear 
approximation Eq.(A2), we have [41]: 


HiT) — -{e\V—SM ^5‘le) — ^ 

r 


(A4) 

with Q = 1^ — S{S*S + ^5*, and = e-y for 

r = {7 o (f,i)}. 


APPENDIX 


A. Formalism of elastic energy 


The energy H{T) of a given spring configuration T = 
{7 O (i, j)} is defined in Eq.(l) as a minimization on the 
positions of the nodes. This minimum can be calculated 
using conjugate gradient methods. However, for small 
mismatches e, it is more efficient to use linear algebra [41], 
as we now recall. Consider a displacement field SRi = 
Ri — Rio, where Rio is the position of the node i in the 
crystal described in the previous section. We define the 
distance [[Rio ~ ^ioll = ''’(i.j)- ^.t first order in SRi, the 
distance among neighboring nodes can be written as: 

\\Ri — RjW = +^^S(^ij'jkSRk + o{SlP) (Al) 

k 

Where S is the structure matrix, which gives the linear 
relation between displacements and changes of distances. 


B. Density of states 

We have shown the density of states converges to 
the one of mean-field networks [42]. Cooling strongly 
suppresses low-frequency vibrational modes, as seen in 
Fig. Al. This temperature effect on the density of states 
is primarily induced by the weak interactions: the den¬ 
sity of states changes little under cooling when a = 0 , 
as appeared in comparing (a) and (b) of Fig. A2. The 
slight change indicates that the density of states depends 
on the presence of redundant constraints. However, when 
a > 0 , the low-temperature density of states strongly dif¬ 
fers from its high-temperature counterpart, as shown in 
Fig. A2(a) and (c). 

The modes that rarefy under cooling are local¬ 
ized vibrations. The participation ratio, P{uj) = 
l^(Ei J2i quantifies the extensity of charac¬ 

teristic modes: P ^ 0 corresponds to a localized mode, 
while P ^ 1 means that the mode extends over the sys¬ 
tem. Both the low and high-frequency ends of the den¬ 
sity of states are reduced under cooling, but the modes 
in the middle are enhanced, as shown in the right panel 
of Fig. Al. This agrees with the small participation ratio 
of modes with low and high frequencies, see Fig. A2(d). 
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In fact, all modes become extended - the participation 
ratio increases over the whole spectrum - when the tem¬ 
perature decreases, as shown in the inset of Fig. Al. 

In addition to localization, another prominent feature 
of reduced low-frequency modes is the power-law diverg¬ 
ing density of states D{uj) ^ see Fig. A2. The 

abundance of low-frequency localized modes appearing 
with a power law density of states signals the “fractons” 
that appear near the rigidity percolation [47, 65, 66]. The 
exponent of the diverging tail, in Fig. A2(a), implies the 
fracton dimension d ~ 0.75, which is consistent with 0.78 
observed for the rigidity percolation [66, 67]. Different 
fracton dimensions d are observed for different coordi¬ 
nation numbers in the case of rigidity window shown in 
Fig. A2(b), although more work would be needed to es¬ 
tablish this fact empirically. 

We discuss when the temperature affects the mode 
with frequency uj in Appendix Sec. C and show illus¬ 
trations of “fractons” in Appendix Sec. D. 


C. Adaptation effects on density of states 

When a > 0, following Eq.(6), we find out the typi¬ 
cal elastic energy corresponding to a mode of frequency 
u! scales as ajiitp' -fa), which is proportional to a for 
w ^ 1, while proportional to 1 when w <C This 

implies that the elastic energy in the degrees of freedom 
corresponding to the modes of low-frequency is of the 
same magnitude as the one in the redundant constraints. 
Similar to the redundant constraints, these low-frequency 
modes are reduced under cooling. 

From Eq.(12b), T*(a;,a) ~ a/(a;^ -f a) gives an esti¬ 
mate of the temperature scale the mode w begins to be 
reduced. The adaptation effect at this temperature scale 
can be seen in the right panel of Fig. Al. For example. 




FIG. Al: (Color online) Changes of density of states 
T) with temperature for the same 2 = —0.055, 
a = 0.0003. Left: density of states in log-log scale. Right: 
density of states normalized by its T = oo value, 
emphasizing its difference under cooling. Inset: participation 
ratio P(a;, T) variation under cooling. 


the dashed green line at T 0.04 <C 1 shows a den¬ 
sity of states with frequencies ui < « 0.01 strongly 

suppressed, while the shape of the density of states with 
w « 0.1 and above is almost unchanged. The dotted pur¬ 
ple line, T Ri 10“^ ~ a, shows a density of states whose 
highest frequency a; ~ 1 is also significantly reduced. 


D. Fractons 


“Fractons” are different from either the low-frequency 
Debye modes or the anomalous modes on the boson peak, 
as shown in Fig. A3. They (Fig. A3(c)) are localized 
and random compared to the Debye modes (Fig. A3 (a)), 
and concentrated on the fractal sets with sharp bound¬ 
aries, unlike the extended anomalous modes (Fig. A3(b)). 
The “fractons” are associated with the collective motion 
of large isostatic or nearly isostatic regions as shown in 
Fig. A4. 


E. Vibrational entropy contribution 


The structure the elastic potential evolve with temper¬ 
ature in the liquid phase of the adaptive network model. 
Freezing into a glass phase eliminates this variability and 
leads to a contribution to the jump of specific heat [68]. 
Our model currently ignores the vibrational part of the 
specific heat, which incorporates that the shape of the 
inherent structure evolves with temperature - not only 
its bottom energy. We estimate this contribution from 
vibrations in this subsection and argue that is is not sig¬ 
nificant for the models we consider. 

The vibrational entropy includes both linear w > 0 and 
floppy w = 0 vibration modes [68]: 

/■ eT 

Svib{T) = [1 - nr{T)] / da;D(a;, ^ + /(’^) 

(A5) 

A sets a cutoff volume for floppy modes, which is approx¬ 
imately the atomic spacing measured in the Lindemann’s 
length: A « (1/0.15)^^ [58], of order 10^ in 3D [57]. / is 
the floppy mode density, dual to the number density of 
redundant constraints flT) = —dz/z nr{T) and thus 
df{T)/dT = dririT) / dT . The jump of specific heat fol¬ 
lows: 
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FIG. A2: (Color online) Density of states D{uj,T) for adaptive networks with different (a) Random diluted networks 
T = oo; a power law D{u}) ~ is shown in the low-frequency range of networks near Zcen- (b) Adaptive networks 

without weak constraints (a = 0) at T = 0.0003; power laws with different exponents are shown for networks in the rigidity 
window: D(uj) ~ for 5z = —0.055, D{uj) ~ ® for 5z = 0.0. (c) Adaptive networks with weak constraints 

(a = 0.0003) at T « a; away from isostatic, the densities of states are gapped between zero frequency and Boson peak, where 
D{uj) ~ u)°. Inset (d) is the participation ratio P{uj,T) at T = oo, see text for definition. 



FIG. A3: (Color online) Vector plots of vibrational modes in randomly diluted networks, N = 100 x 100. (a) A typical Debye 
mode, Sz = 0.501, uj = 0.017. (b) A typical anomalous mode on boson peak, 5z = —0.049, oj — 0.011. (c) A typical fracton, 

Sz = -0.049, LJ = 0.0007. 




^ dnr{T) 
9 QJ. 


In A — 


dujDTg (oj) In 


The derivatives on InT in Eq.(A5), continuous at the 
glass transition, have been subtracted. 

We estimate the upper limit of the vibrational contri¬ 
bution. (1) The first term in Eq.(A6): Debye frequency 
(jJd sets the upper limit of the integral in the bracket, 
— ln(eTg/fiuJD)- As the glass transition temperature Tg 
and Debye temperature du = huJo/ks are usually of 
the same order, the bracket in the first term is doini- 


hu! 


+ [1 - nr{Tg)] / dw Tj 


dDriuj) 


dT 


In^ 

OLO 


(A6) 


nated by InA. From Eqs.(ll), we have dnr/dhiT\Tg ~ 
^nex(Tg) < ^noy/T^ < 0.02i/a, and InA fa 5 in 2D. 
Compared to the specific heat values, which are of or¬ 
der one shown in Fig. 7, and the scalings of the minima 
—O.l/liiQ; given in [41], the contribution, is in¬ 

significant if 0 < a < 0.1. 

(2) The second term in Eq.(A6): The upper limit of 
the bracket is 1. Replacing ln(eT/ficj) with its upper 
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limit In A, we simplify the integral to f dujTdD/dT. We 
can estimate the upper limit of the derivative in the 
integral approximately by Aut / A In T, where Any is 
the number density of the modes reduced under cool¬ 
ing. Aht ~ 0.2 « 0.01, roughly the num¬ 

ber fraction of “fractons” suppressed under cooling. To¬ 
gether, the upper limit of the contribution of the second 
term is An^/ In 10 x In A Ri 0.03, which is moderate com¬ 
pared to the values of order one. 

Therefore, the vibrational entropy contributes mildly 
to the jump of specific heat and does not change the qual¬ 
itative behavior of Ac in our model of network glasses. 
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